Data Comparison - NES vs SES (WIP)

0.1 Import Data

# load wealingNES package
library(weanlingNES)

# load dataset
data("data_nes")
data("data_ses")

# merge into one dataset
data_comp <- rbind(
  rbindlist(data_nes$year_2018) %>%
    .[, sp := "nes"],
  rbindlist(data_ses$year_2014) %>%
    .[, sp := "ses"],
  use.names = T,
  fill = T
)

0.2 Plots

Let’s try to recreate the same plots Grecian, et al. (2022) did to compare the evolution of some behavioral parameters across time between nes and ses.

0.2.1 Max Depth

p1 <- plot_comp(data_comp[phase == "day", ], "maxdepth", nb_days = 300) +
  labs(
    x = "# days since departure",
    y = "Maximum Depth (m)",
    title = "Day"
  ) +
  scale_y_continuous(limits = c(0, data_comp[, max(maxdepth)])) +
  theme_jjo()
p2 <- plot_comp(data_comp[phase == "night", ], "maxdepth", nb_days = 300) +
  labs(
    x = "# days since departure",
    y = "Maximum Depth (m)",
    title = "Night"
  ) +
  scale_y_continuous(limits = c(0, data_comp[, max(maxdepth)])) +
  theme_jjo() +
  theme(
    axis.title.y = element_blank(),
    axis.ticks.y = element_blank(),
    axis.text.y = element_blank()
  )
ggarrange(p1, p2, ncol = 2, common.legend = TRUE, legend = "bottom")

Figure 1: Estimated temporal changes in maximum depth (m)

0.2.2 Dive Duration

p1 <- plot_comp(data_comp[phase == "day", ], "dduration", nb_days = 300) +
  labs(
    x = "# days since departure",
    y = "Dive Duration (s)",
    title = "Day"
  ) +
  scale_y_continuous(limits = c(0, data_comp[, quantile(dduration,0.99)])) +
  theme_jjo()
p2 <- plot_comp(data_comp[phase == "night", ], "dduration", nb_days = 200) +
  labs(
    x = "# days since departure",
    y = "Dive Duration (s)",
    title = "Night"
  ) +
  scale_y_continuous(limits = c(0, data_comp[, quantile(dduration,0.99)])) +
  theme_jjo() +
  theme(
    axis.title.y = element_blank(),
    axis.ticks.y = element_blank(),
    axis.text.y = element_blank()
  )
ggarrange(p1, p2, ncol = 2, common.legend = TRUE, legend = "bottom")

Figure 2: Estimated temporal changes in dive duration (s)

0.2.3 “bADL”

# data nes
days_to_keep_nes = data_comp[sp=="nes",
                                .(nb_dives = .N),
                                by = .(.id, day_departure)] %>%
  .[nb_dives >= 50,]
# keep only those days
data_nes_complete_day = merge(data_comp[sp == "nes",],
                                      days_to_keep_nes,
                                      by = c(".id", "day_departure"))
# data ses
days_to_keep_ses = data_comp[sp=="ses",
                                .(nb_dives = .N),
                                by = .(.id, day_departure)] %>%
  .[nb_dives >= 8,]
# keep only those days
data_ses_complete_day = merge(data_comp[sp == "ses",],
                                      days_to_keep_ses,
                                      by = c(".id", "day_departure"))
# data plot
dataPlot = rbind(data_nes_complete_day[divetype=="1: foraging",
                                         .(badl = quantile(dduration, 0.95)),
                                         by = .(.id, day_departure, sp)],
                 data_ses_complete_day[divetype=="1: foraging",
                                         .(badl = quantile(dduration, 0.95)),
                                         by = .(.id, day_departure, sp)])

# comparative plot
plot_comp(dataPlot, "badl", nb_days = 300, alpha_point = .1) +
  labs(
    x = "# days since departure",
    y = "bADL (s)",
    title = "Day"
  ) +
  scale_y_continuous(limits = c(0, dataPlot[, quantile(badl,0.99)])) +
  theme_jjo()

Figure 3: Estimated temporal changes in bADL (s)

0.2.4 Drift Rate

# calculate drift rate per day, id and sp
dataPlot = data_comp[divetype == "2: drift",
  .(driftrate = quantile(driftrate, 0.5)),
  by = .(day_departure, .id, sp)
]

# comparative plot
p1 = plot_comp(dataPlot, "driftrate", nb_days = 300, alpha_point = .1) +
  labs(
    x = "# days since departure",
    y = "Drift Rate (m/s)",
    title = "Day"
  ) +
  theme_jjo() +
  theme(legend.position = "bottom")
p2 = ggplot(dataPlot, aes(x = day_departure, y = driftrate, col = .id)) +
  geom_point(show.legend = FALSE) +
  facet_grid(sp~.) +
  theme_jjo() +
  theme(axis.title.y = element_blank())
ggarrange(p1, p2, ncol = 2)

Figure 4: Estimated temporal changes in drift rate (m/s)

It’s weird that there is still a bimodality in the maxdepth’s distribution for northern elephant seal, even after splitting day and night.

To investigate why is that, let’s try several representation of this data for the individual 2018070 (nes).

# let's pick an individual
data_test <- data_comp[.id == "ind_2018070", ]

# first we average `lightatsurf` by individuals, day since departure and hour
dataPlot <- data_test[, .(lightatsurf = median(lightatsurf, na.rm = T),
                          phase = first(phase)),
                      by = .(.id,
                             day_departure,
                             date = as.Date(date),
                             hour)]

# then we choose the variable to represent
i <- "maxdepth"
ggplot(data = melt(data_test[, .(.id, date, get(i), phase)],
                   id.vars = c(".id",
                               "date",
                               "phase")),
       aes(x = as.Date(date),
           y = value,
           col = phase)) +
  geom_point(alpha = 1 / 10,
             size = .5) +
  facet_wrap(. ~ .id, scales = "free") +
  scale_x_date(date_labels = "%m/%Y") +
  labs(x = "Date", y = i) +
  theme_jjo() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "bottom") +
  guides(colour = guide_legend(override.aes = list(size = 7,
                                                   alpha = 1)))

Figure 5: Evolution of the maximum depth reached across time for the individual 2018070

For this individual we can see that there are still a lot of shallow dives during the day. I first thought it was because this individual would probably have started reaching shallow water earlier in the day, at dusk, rather than waiting for complete darkness. To test this hypothesis, I tried to visualize the maxdepth over an entire day throughout the trip.

# this is art!
htmltools::tagList(list(
plot_ly() %>%
  add_trace(
    data = data_test,
    x = ~hour,
    y = ~day_departure,
    z = ~ -maxdepth,
    marker = list(
      size = 1,
      opacity = 0.5
    ),
    mode = "markers",
    type = "scatter3d",
    text = ~divenumber
  ) %>%
  add_trace(
    x = ~hour,
    y = ~day_departure,
    z = ~ (lightatsurf / 200) * 20,
    intensity = ~ lightatsurf / 200,
    data = dataPlot,
    type = "mesh3d",
    showlegend = FALSE
  ) %>%
  layout(
    scene = list(
      zaxis = list(title = "Maximum depth (m)"),
      yaxis = list(title = "# days since departure"),
      xaxis = list(title = "Hour")
    ),
    legend = list(itemsizing = "constant")
  )))

We can see the bimodality of maxdepth within a day is not related to the time of the day, since swallow and deep dives might both occurred at the same time in the day. Let’s see what happen if we color dives with divetype.

# this is art!
htmltools::tagList(list(
  plot_ly() %>%
  add_trace(
    data = data_test,
    x = ~hour,
    y = ~day_departure,
    z = ~ -maxdepth,
    color = ~divetype,
    marker = list(
      size = 1,
      opacity = 0.5
    ),
    mode = "markers",
    type = "scatter3d",
    text = ~divenumber
  ) %>%
  add_trace(
    x = ~hour,
    y = ~day_departure,
    z = ~ (lightatsurf / 200) * 20,
    intensity = ~ lightatsurf / 200,
    data = dataPlot,
    type = "mesh3d",
    showlegend = FALSE
  ) %>%
  layout(
    scene = list(
      zaxis = list(title = "Maximum depth (m)"),
      yaxis = list(title = "# days since departure"),
      xaxis = list(title = "Hour")
    ),
    legend = list(itemsizing = "constant")
  )
))

Well it seems that the bimodality observed in the distribution of maxdepth is mostly explained by divetype, with transit dives occurring at deeper dives that foraging and drift dives (at least for northern elephant seal). We can actually look at the same result with a simple 2D plot:

ggplot(
  data = melt(data_test[, .(.id, date, get(i), divetype)],
              id.vars = c(
                ".id",
                "date",
                "divetype"
              )
  ),
  aes(
    x = as.Date(date),
    y = value,
    col = divetype
  )
) +
  geom_point(
    alpha = 1 / 10,
    size = .5
  ) +
  facet_wrap(. ~ .id, scales = "free") +
  scale_x_date(date_labels = "%m/%Y") +
  labs(x = "Date", y = i) +
  theme_jjo() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.position = "bottom"
  ) +
  guides(colour = guide_legend(override.aes = list(
    size = 7,
    alpha = 1
  )))

Figure 6: Kind of the same graph, but without looking at the hour of the day

And in bonus, we can add the bathymetry to the 3D plot.

htmltools::tagList(list(
  plot_ly() %>%
  add_trace(
    data = data_test,
    x = ~hour,
    y = ~day_departure,
    z = ~ -maxdepth,
    color = ~divetype,
    marker = list(
      size = 1,
      opacity = 0.5
    ),
    mode = "markers",
    type = "scatter3d",
    text = ~divenumber
  ) %>%
  add_trace(
    x = ~hour,
    y = ~day_departure,
    z = ~ (lightatsurf / 200) * 20,
    intensity = ~ lightatsurf / 200,
    data = dataPlot,
    type = "mesh3d",
    showlegend = FALSE
  ) %>%
  add_trace(
    x = ~hour,
    y = ~day_departure,
    z = ~bathy,
    data = data_test,
    type = "mesh3d"
  ) %>%
  layout(
    scene = list(
      zaxis = list(title = "Maximum depth (m)"),
      yaxis = list(title = "# days since departure"),
      xaxis = list(title = "Hour")
    ),
    legend = list(itemsizing = "constant")
  )))

What about Southern Elephant Seals ?!?

To investigate why is that, let’s try several representation of this data for the individual 140059 (ses).

# let's pick an individual
data_test <- data_comp[.id == "ind_140059",]

# first we average `lightatsurf` by individuals, day since departure and hour
dataPlot <- data_test[, .(lightatsurf = median(lightatsurf, na.rm = T),
                          phase = first(phase)),
                      by = .(.id,
                             day_departure,
                             date = as.Date(date),
                             hour)]

# then we choose the variable to represent
i <- "maxdepth"
ggplot(data = melt(data_test[, .(.id, date, get(i), phase)],
                   id.vars = c(".id",
                               "date",
                               "phase")),
       aes(x = as.Date(date),
           y = value,
           col = phase)) +
  geom_point(alpha = 1 / 5,
             size = .5) +
  facet_wrap(. ~ .id, scales = "free") +
  scale_x_date(date_labels = "%m/%Y") +
  labs(x = "Date", y = i) +
  theme_jjo() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "bottom") +
  guides(colour = guide_legend(override.aes = list(size = 7,
                                                   alpha = 1)))

Figure 7: Evolution of the maximum depth reached across time for the individual 140059

For this individual we can see that until the half of its trip, there is no distinction between day and night in terms of maxdepth. But After, he starts to dive deeper during the day, and shallower during the night. Contrary to the previous nes individual, there seems to be a clear pattern. Let’s look at the graph in 3D:

# this is art!
htmltools::tagList(list(
  plot_ly() %>%
  add_trace(
    data = data_test,
    x = ~hour,
    y = ~day_departure,
    z = ~ -maxdepth,
    marker = list(
      size = 1,
      opacity = 0.8
    ),
    mode = "markers",
    type = "scatter3d",
    text = ~divenumber
  ) %>%
  add_trace(
    x = ~hour,
    y = ~day_departure,
    z = ~ (hour * 0.1) + 20,
    intensity = ~phase_bool,
    data = dataPlot[][, phase_bool := fifelse(phase == "night", 0, 1)][],
    type = "mesh3d",
    showlegend = FALSE
  ) %>%
  layout(
    scene = list(
      zaxis = list(title = "Maximum depth (m)"),
      yaxis = list(title = "# days since departure"),
      xaxis = list(title = "Hour")
    ),
    legend = list(itemsizing = "constant")
  )))

We can clearly see that pass half of its trip, this individual starts to dive deeper during the day. Let’s see what happen if we color dives with divetype.

# this is art!
htmltools::tagList(list(
  plot_ly() %>%
  add_trace(
    data = data_test,
    x = ~hour,
    y = ~day_departure,
    z = ~ -maxdepth,
    color = ~divetype,
    marker = list(
      size = 1,
      opacity = 0.8
    ),
    mode = "markers",
    type = "scatter3d",
    text = ~divenumber
  ) %>%
  add_trace(
    x = ~hour,
    y = ~day_departure,
    z = ~ (hour * 0.1) + 20,
    intensity = ~phase_bool,
    data = dataPlot[][, phase_bool := fifelse(phase == "night", 0, 1)][],
    type = "mesh3d",
    showlegend = FALSE
  ) %>%
  layout(
    scene = list(
      zaxis = list(title = "Maximum depth (m)"),
      yaxis = list(title = "# days since departure"),
      xaxis = list(title = "Hour")
    ),
    legend = list(itemsizing = "constant")
  )
))

It’s hard to tell something, most of the dives seem to be foraging dives, and follow the pattern identified previously. The other interesting results would be a lot of benthic dives at the beginning, and drift dives occurring mostly during the day. Let’s look at a simple 2D plot:

ggplot(
  data = melt(data_test[, .(.id, date, get(i), divetype)],
              id.vars = c(
                ".id",
                "date",
                "divetype"
              )
  ),
  aes(
    x = as.Date(date),
    y = value,
    col = divetype
  )
) +
  geom_point(
    alpha = 1 / 5,
    size = .5
  ) +
  facet_wrap(. ~ .id, scales = "free") +
  scale_x_date(date_labels = "%m/%Y") +
  labs(x = "Date", y = i) +
  theme_jjo() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.position = "bottom"
  ) +
  guides(colour = guide_legend(override.aes = list(
    size = 7,
    alpha = 1
  )))

Figure 8: Kind of the same graph, but without looking at the hour of the day

And in bonus, we can add the bathymetry to the 3D plot.

htmltools::tagList(list(plot_ly() %>%
  add_trace(
    data = data_test,
    x = ~hour,
    y = ~day_departure,
    z = ~ -maxdepth,
    color = ~divetype,
    marker = list(
      size = 1,
      opacity = 0.8
    ),
    mode = "markers",
    type = "scatter3d",
    text = ~divenumber
  ) %>%
  add_trace(
    x = ~hour,
    y = ~day_departure,
    z = ~ (lightatsurf / 200) * 20,
    intensity = ~ lightatsurf / 200,
    data = dataPlot,
    type = "mesh3d",
    showlegend = FALSE
  ) %>%
  add_trace(
    x = ~hour,
    y = ~day_departure,
    z = ~bathy,
    data = data_test,
    type = "mesh3d"
  ) %>%
  layout(
    scene = list(
      zaxis = list(title = "Maximum depth (m)"),
      yaxis = list(title = "# days since departure"),
      xaxis = list(title = "Hour")
    ),
    legend = list(itemsizing = "constant")
  )))